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Abstract 

We constrain the light CP-odd Higgs A® in NMSSM via the rare decay 7r° — > e + e _ . It 
is shown that the possible 3a discrepancy between theoretical predictions and the recent 
KTeV measurement of B(ir° — > e + e~) cannot be resolved when the constraints from 
T — ► a„ and 7r° — ► 77 are combined. Furthermore, the combined constraints also 

exclude the scenario involving m A o = 214.3 MeV, which is invoked to explain the anomaly 
in the S + — > pfi + [i~ decay found by the HyperCP Collaboration. 

PACS Numbers: 13.25.Cq, 14.80.Cp 
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1 Introduction 



Theoretically, the rare decay 7r° — > e + e~ starts at the one loop level in the standard model (SM), 
which has been extensively studied [II [21 El III El El El El El E3II| since the first investigation in QED 
by Drell [T]. It is nontrivial to make precise predictions of the branching ratio Bsm(^° —* e + e~) 
because its sub-process involves the n° — > 7*7* transition form factor. In Refs. [21 El IU E] , the 
decay was studied via the Vector-Meson Dominance (VMD) approach, where the results are in 
good agreement with each other and converge in B(tt° — > e + e~) ~ 6.2 — 6.4 x 1CT 8 . By using 
the measured value of B(r] — > to fix the counterterms of the chiral amplitude in Chiral 

Perturbation Theory (ChPT), Savage et al. predicted B(tt° -> e+e") = (7 ± 1) x 10~ 8 [6]. 
Using a procedure similar to that used in Ref. [H] (although with an updated measurement of 
B{r\ — > Dumm and Pich predicted (8.3±0.4) x 1CT 8 [7\. Alternatively, using the lowest 

meson dominance (LMD) approximation to the large- N c spectrum of vector meson resonances 
to fix the counterterms, Knecht et al. predicted (6.2 ± 0.3) x 10~ 8 [8], which is about 4<r lower 
than the value predicted by Ref. [7] but which agrees with the others. Most recently, using 
a dispersive approach to the amplitude and the experimental results of the CELLO [TT] and 
CLEO [12] Collaborations for the pion transition form factor, Dorokhov and Ivanov [9] have 
found that 

Bsm(tt° -> e + e") = (6.23 ± 0.09) x 10~ 8 , (1) 

which is consistent with most theoretical predictions of Bsm(^° ~* e + e~) in the literature. 
Moreover, their prediction that B(r] — > /i + yU~) = (5.11±0.2) x 10~ 6 agrees with the experimental 
data (which gives a value of (5.8 ± 0.8) x 10" 6 p3]). 

Experimentally, the accuracy of the measurements of the decay has increased significantly 
since the first tt° — > e + e~ evidence was observed by the Geneva-Saclay group [Hj in 1978 with 
Bsm(^° —* e + e~) = (221^) x 10^ 8 . A detailed summary of the experimental situation can be 
found in Ref. [IS]. Recently, using the complete data set from KTeV E799-II at Fermilab, the 
KTeV Collaboration has made a precise measurement of the 7r° — > e + e~ branching ratio [Tfi"] 

B™-™ d (jP -> e + e~) = (7.48 ± 0.29 ± 0.25) x 10~ 8 , (2) 

after extrapolating the full radiative tail beyond (m e + e - /m^o) 2 > 0.95 and scaling their result 
back up by the overall radiative correction of 3.4%. 
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As was already noted in Ref. [5], the SM prediction given in Eq.flT]) is 3.3cr lower than the 
KTeV data. The authors have also compared their result with estimations made by various 
approaches in the literature and found good agreements. Further analyses have found that 
QED radiative contributions [T7] and mass corrections [TS] are at the level of a few percent 
and are therefore unable to reduce the discrepancy. Although the discrepancy might be due to 
hadronic dynamics that are as of yet unknown, it is equally possible that this discrepancy is 
caused by the effects of new physics (NP). In this Letter we will study the latter possibility. 

As is known that leptonic decays of pseudoscalar mesons are sensitive to pseudoscalar weak 
interactions beyond the SM. Precise measurements and calculations of these decays will offer 
sensitive probes for NP effects at the low energy scale. Of particular interest to us is the rare 
decay ir° — > e + e~, which could proceed at tree level via a flavor- conserving process induced by 
a light pseudoscalar Higgs boson A\ in the next-to-minimal supersymmetric standard model 
(NMSSM) [19]. We will look for a region of the parameter space of NMSSM that could resolve 
the aforementioned discrepancy of B(tt° — > e + e~) at la. Then, we combine constraints from 

and the recent searches for T(1S), (3S) — > ^A\ by CLEO [20] and BaBar [2T], respectively. 

2 The amplitude of tt° e+e~ in the SM and the NMSSM 

The NMSSM has generated considerable interest in the literature, which extends the minimal 
supersymmetric SM (MSSM) by introducing a new Higgs singlet chiral superfield S to solve 
the known \x problem in MSSM. The superpotential in the model is [TU] 

W NMSSM = QH u h u U c + H d Qh d D c + H d Lh e E c + XSH u H d + ^S 3 , (3) 

where k is a dimensionless constant and measures the size of Peccei-Quinn (PQ) symmetry 
breaking. 

In addition to the two charged Higgs bosons, H , the physical NMSSM Higgs sector consists 
of three scalars h°, Hf 2 and two pseudoscalars A\ 2 . As in the MSSM, tan (3 = v u j v d is the ratio 
of the Higgs doublet vacuum expectation values v u = (H®) = v sin j3 and v d = (H®) = v cos /?, 
where v = v\ + v\ = \/2mw / ' 9 — l7AGeV. Generally, the masses and singlet contents of the 
physical fields depend strongly on the parameters of the model (such as, in particular, how well 
the PQ symmetry is broken). If the PQ symmetry is slightly broken, then A\ can be rather 
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Figure 1: Relevant Feynman diagram within NMSSM. 



light, and its mass is given by 

mj, = 3™* + O(J^) (4) 

with the vacuum expectation value of the singlet x = (S); meanwhile, another pseudoscalar A® 
has a mass of order of mjj±. 

For 7T° — > e + e~ decay, the NMSSM contributions are dominated by A®. The couplings of 
A\ to fermions are [22] 

ff = [x d m d d^d + X u m u uj 5 u + X^mei^lj A\ (5) 

where X d = Xg = ^5- and X u = X d j tan 2 /3; thus, the contribution of the u^§uA\ term in 
7T° — > e + e~ could be neglected in the large tan (3 approximation. 

To the leading order, the relevant Feynman diagram within NMSSM is shown in Fig. [TJ We 
obtain its amplitude as 

G 1 

M A o = — |m e my^.-j ZwT X d> ( 6 ) 

which is independent of m d , since m d in the coupling of A^dj^d is canceled by the m d term of 
the hadronic matrix 

<0|W>^/„.|fc. (7) 
In the SM, the normalized branching ratio of ir° — > e + e~ is given by [§] 

«<«° - « + - %^ = 2 (?^) 2 ^K.)M«.)l 2 (8) 



where P e (m? a ) = J\ — 4^r- and A(m? ) is the reduced amplitude. 

To add the NMSSM amplitude to the above amplitudes consistently, we rederive the SM 
amplitude to look into possible differences between the conventions used in our Letter and the 
ones used in Ref. [2j • The Feynman diagram that proceeds via two photon intermediate states 
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Figure 2: Triangle diagram for 7r° — ► e + e process. 



is shown in Fig. [21 We start with the 7T°7*7* vertex 



= -i e 2 e^p k a (q - k) 13 f rr F n o rr (k 2 , (q - kf 



(9) 



where k and q — k are the momenta of the two photons, / 7 * 7 * = 47r2 is the coupling constant 
of 7r° to two real photons. F n o 1 * 1 *(k 2 , (q — k) 2 ) is the transition form factor 7T° — > 7*7*, which 
is normalized to F 7r o 7 * 7 *(0, 0) = 1. The amplitude of Fig. [2] is written as 

o /* d k 



/;;/ 



(2vr) 4 (A; 2 + ie) ((fc - g) 2 + j e ) ((fc - p ) 2 - m e + ie) 



(10) 



with 



IT = u(p,s)Y(]t- ft +m e )Yv(q-p,s'). (11) 

There is a known, convenient way to calculate L^ u with the projection operator for the outgoing 
e + e~ pair system [23] 



V(q-p,p) 



V2 



[v(q — p, +) ® u(p, -) + u(g - p, -) (g) +)] 
^= [ - 2rJ W yV + |wr(^(9 " P) T - (g - + rf] (12) 



where t = q 2 = m 2 . After some calculations, we get 



Msm(k — > e + e ) = 2v / 2a 2 m e m 7r o / 7 * 7 * A(m^) 



(13) 



where the reduced amplitude A(g 2 ) is 



. , 9. 2 z f d k 

Aq 2 ) = -o / -r 



/c 2 g 2 - (g ■ fc)' 



7r 2 (A; 2 + ie) ((fc — g) 2 + ze) ((fc — p) 2 — m e + ie) 



^F^o rr {k 2 ,{q-k) 2 ). (14) 



We note that the -A(g 2 ) derived here is in agreement with Ref. pj. Further evaluation of the 
integrals of «4(g 2 ) is quite subtle and lengthy [21 [21], and only the imaginary part of A(m 2 ) 



can be obtained model- independently [TJ [2]. In the following calculations, we quote the result 
of Ref. [9], 

A{ml) = (10.0 ± 0.3) - il7.5. (15) 
With Eq. (EJ) and Eq. ( fTBl . we get the total amplitude 



M. = 2v / 2a 2 m e m T o / 7 * 7 * ^4(^) T=Tn e m^o f^o 



s/2 



(16) 



3 Numerical analysis and discussion 

Now, we are ready to discuss the effects of A® numerically, with a focus on the m A o < 2mt, 
scenarios. The dependence of B(ir° — > e + e~) on the parameter \Xd\ is shown in Fig. [3] with 
m A o = m n /2, 214.3MeV, 3GeV as benchmarks. We have used the input parameters B(tt° — > 
77 ) = 0.988 and f n o = (130.7 ± 0.4) MeV [33]. As shown in Fig. El B(7i° -> e + e") is 
very sensitive to the parameter \X^\ and m A o. For m A o < m n o, the NMSSM contribution is 
deconstructive and reduces B(ir° — > e + e~) at small \Xj\ region. For m A o > m n o, the NMSSM 



contribution is constructive and could enhance B(n c 



e T e" 



to be consistent with the KTeV 



measurement S^ o Te ™ d (7r 



m A o). 



e + e 



(7.48 ± 0.38) x 10 8 (where \X d \ strongly depends on 




Figure 3: The dependence of B(ir° —>■ e + e~) on the parameter \Xd\ with m A o = m n o/2, 
214.3 MeV and 3 GeV, respectively. The horizontal lines are the KTeV data, where the solid 
line is the central value and the dashed ones are the error bars (lex). 

I. Constraint on the scenario of m A o = 214.3MeV 

It is interesting to note that the HyperCP Collaboration [25] has observed three events for 
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the decay E + — » pn + fj. with a narrow range of dimuon masses. This may indicate that the 
decay proceeds via a neutral intermediate state, S + — > pP°,P° — > with a P° mass of 

214.3 ± 0.5MeV. The possibility of P° has been explored in the literature [26! EEJ EHJ [30] . The 
authors have proposed A\ as a candidate for the P°, and have also shown that their explanation 
could be consistent with the constraints provided by K and B meson decays [26l [27] . It would 
be worthwhile to check on whether the explanation could be consistent with the tt° — > e + e~ 
decay. 

Taking m A o = 214.3 MeV, we find that B(tt° — > e + e~) is enhanced rapidly and could be 
consistent with the KTeV data within la for 

\X d \ = 14.0 ± 2.4. (17) 

However, the upper bound \X d \ < 1.2 from the a M constraint has been derived and used in the 
calculations of Ref. [2"6"l |2"9"] . So, with the assumption that m A o = 214.3 MeV, our result of \Xd\ 
violates the upper bound with a significance of 5a. 

Recently, CLEO [20] and BaBar [21] have searched for the CP-odd Higgs boson in radiative 
decays of T(IS') — > and T(3S) — > 'yA®, respectively. For m A o = 214MeV, CLEO gives the 
upper limit 

B(T(1S) -> iA\) < 2.3 x 10" 6 (90% C.L.) (18) 

which constrains \X&\ < 0.16. 

The BaBar Collaboration has searched for A\ through T(3S*) — > 'yA®, A\ — * invisible in the 
mass range m A o < 7.8 GeV ^T\. From Fig. 5 of Ref. [21], we read 

B(T(3S) -> 7 A?) x B{A\ -> invisible) < 3.5 x 10~ 6 (90% C.L.) (19) 

for m A o = 214 MeV. Assuming B(A\ — > invisible) ~ 1, we get the conservative upper limit 
\X d \ < 0.19 

All of these upper limits are much lower than the limit of Eq{T7Jset by 7r° — > e + e~; therefore, 
the scenario where m A o ~ 214 MeV in NMSSM could be excluded by combining the constraints 
from 7T° — > e + e~ and the direct searches for T radiative decays. 
II. Constraints on the parameter space of m A o — \X d \ 

To show the constraints on NMSSM parameter space from n° — > e + e~, we present a scan of 
m A o — \Xd\ space, as shown in Fig. [H In order to scan the region of m A o ~ m^o, the amplitude 
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Figure 4: Constraints on the NMSSM parameter space through B(ir°^e + e ), B(T(1S) — > jA®), 
B(Y(3S) — > 7^5) and a M respectively. The shaded regions are allowed by the labeled processes. 



of the A\ contribution in Eq. ([6]) is replaced by the Breit-Wigner formula 



M a°, = - -7=m c m,o/x»- 2 — 

1 V2 m Z o - m 2 



(20) 



o A v +iY(A\)m A o 

With the assumption that A\ just decays to electron and photon pairs for m A o ~ m n o, the 
decay width of A\ could be written as 



rK) = r(yi;^eV) + r(A^ 77 ) 



(21) 



with 



e e 



r(A? - 77) 



•F 2 v-2 /-i m e 



2 

u i ' 



(22) 



where r = 1 for leptons and r = N c for quarks, k\ = m^/m A o and Qi is the charge of the 
fermion in the loop. The loop function F{kj) reads 

i \ 2 

F(ki 



for fct > j, 



for < 4 • 



As shown in Fig. HI only two narrow connected bands of the \Xd\ — space survive after 
the KTeV measurement of B(7r°^e + e~), which show that 7r°^e + e~ is very sensitive to NP 
scenarios with a light pseudoscalar neutral boson. 

In the following, we will determine which part of the remaining parameter space could 
satisfy the constraints enforced by radiative T decays and a M simultaneously. 



S 



To include the a M constraint, we use the experimental result that [32] a^Exp) = (1 1659208. 0± 
6.3) x KT 10 and the SM prediction [33] a^SM) = (11659177.8 ±6.1) x 1(T 10 . The discrepancy 
is 

A a M = a^Exp) - a^SM) = (30.2 ± 8.8) x 10 _10 (3.4o-) (23) 

which is established at a 3.4a level of significance. 
The contributions of A\ to a M are given by [31] 



with 




It has been found that the A\ contribution is always negative at the one loop level and worsens 
the discrepancy in a M ; however, it could be positive and dominated by the two loop contribution 
for A± > 3GeV[3l]. One should note that there are other contributions to a M in NMSSM; for 
instance, the chargino/sneutino and neutralino/smuon loops. Moreover, the discrepancy Aa^ 
could be resolved without pseudoscalars [M]. So, putting a constraint on \Xd\ via a M is a rather 
model-dependent process. There are two approximations with different emphases on the role 
of namely, (i) assuming that Aa M is resolved by other contributions and requiring that 
A\ contributions are smaller than the la error-bar of the experimental measurement, and (ii) 
assuming that the A\ contributions are solely responsible for Aa M . In Ref. [26], approximation 
(i) has been used to derive an upper bound of \X d \ < 1.2. We present the a M constraints with 
the two approximations which are shown in Figs. H^a) and (b), respectively. 

From Fig. H(a), we can find that there are two narrow overlaps between the constraints 
provided by a M and B(ir°— >e + e~): one is for m A o ~ 3 GeV with \Xd\ > 150 and another one is 
for m A o ~ 135 MeV with \X d \ < 1. 
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In the searches for T — > '-/A® decays, CLEO [20] obtains the upper limits for the product of 
B(T(1S) — > 7 ^4?) and B(A\ — > r + r _ ) or B(A\ — > fi + u~), while BaBar presents upper limits on 
B{T(3S) -> 7 A?) x B(A\ -> invisible). All these limits fluctuate with the mass of frequently. 
For simplicity, we take the loosest upper limit B(T(1S) — > 7^4?) x £>(A° — > t + t~) < 6 x 10~ 5 
of CLEO and assume B(A\ — > r + r _ ) = 1. Similarly, we also use the loosest upper limits 
on B(T(3S) -> 7A?) x -> invisible) < 3.1 x 10~ 5 of BaBar [21J and assume B{A\ -> 

invisible) = 1. With the loosest upper limits, we get their bounds on the |JT d | — m^o space, 
which are shown in Fig. HI From the figure, we can see the bounds (excluding the parameter 
space Xd > 1) for < m A a < 7.8 GeV. Fig. IU(b) shows that there is no region of parameter 
space satisfying all the aforementioned constraints if the contribution of A\ is required to solely 
resolve the a M discrepancy. 

Of particular interest, as shown in Fig.H](a), is the parameter space around m A o ~ 135 MeV 
with \Xd\ < 1 (which is still allowed with approximation (i)). To make a thorough investigation 
of the space, we read off the upper limits of BaBar [21] from Fig. 5 for the value m A o ~ 
135 MeV: B(T(3S) -> 7A?) x B{A\ -> invisible) < 3.3 x 10~ 6 . With the assumption that 
B(A® — > invisible) ~ 1 and the constraints from B(tt°— >e + e~), we get 

\X d \ = 0.10 ±0.08, m^o = 134.99 ± 0.01 MeV, (26) 

where the constraint on m A o is dominated by B(n — >e + e~) and the limit of \X^\ is dominated 
by B(T(3S) — > 7^4?). At first sight, the uncertainties in the abovementioned two parameters 
are too different. We find that the difference arises from our assumption T(Ai) ~ T(Ai — > 
e + e~) + T(A® — > 77). From Eqs. fT20l and fl2H) . one can see that the Xj factor in M. A o 
could be canceled out by the one in T(Al) when m A o approaches m n o, which results in a very 
sharp peak for position of m A a. Thus, with the well measured quantities given in Eq. ( 1201 ) 
and the sensitivity of the peak, m A o turns out to be well-constrained. Furthermore, if we take 
m A o = m^o, we find that X\ is canceled out exactly, so there is no parameter to tune; however, 
we have B(ir — > e + e~) ^> 1, which violates the unitary bound and is thus excluded. 

From the results of Eqj26l we obtain 5a^{A\) = (-9.2 ± 8.9) x 10~ 12 with tan/3 = 30 
as a benchmark, which is small enough to be smeared by the chargino/sneutrino and neu- 
tralino/smuon contributions. Moreover, we have 

r(A?) = (5.7 ± 5.5) x 10~ 13 MeV, (27) 
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Figure 5: sin 29 versus the mass difference of the unmixed states with \X^\ = 0.05 and 
0.18. The solid and the dashed lines denote the real and the imaginary parts of sin 29, 
respectively. 

which corresponds to t(A®) ~ 1.2 x 10~ 9 s (cr ~ 36 cm). 

For the case where A\ decays mostly to invisible particles, we take the width of A® as a free 
parameter and get r(A?) < 8.24 x 10~ 6 GeV, m A o = 134.99 ± 0.02 MeV and \X d \ < 0.18. In 
this case, m A o can equal m^o, and it is found that T(A^) < 3.3 x 10~ 6 GeV and \Xd\ < 0.18. 
III. The resonant effects of m A o ~ m n o 

So far we have included only the width effects of A\ with the Breit-Wigner formula for the 
propagator of A\. When the masses of A\ and 7r° are very close, the mixing between the two 
states could modify the parton level 7r° — A\ coupling. In a manner analogous to Ref . [35] , the 
mixing can be described by introducing off-diagonal elements in the A\ — 7r° mass matrix 



M' 



m 



A a - im A oT 
5m 2 



m 2 




(28) 



with 5m 2 = y / G F /4y / 2A 

o??7, nXd. The complex mixing angle 9 between the states is given by 

(5m 2 ) 2 



sin^ 29 



IK? 



1( ) - m 2 - im A oT A o + im^T^o) 2 + (5m 
The mass eigenstates Af and 7r /0 are obtained as 

A'? = ^(A lC os9 + 7i°sm9), 
tt /0 = — (-A°,sm9 + tt°cos9), 



2\2 



(29) 



(30) 
(31) 



where N = ^/| sin 6*| 2 + | cos^l 2 . Then, we can write the decay amplitude of the "physical" 
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state 7T as 

\M(tt'° ^ e + e~)\ 2 = ^(|cos#| 2 |.M(7r ^e + e-)| 2 + | sm^| 2 |>l(A? e + e")| 2 ). (32) 

Obviously, we obtain the SM result when 9 is small. 

With \Xd\ = 0.05 and 0.18, Fig. [5] shows sin26 l as a function of the difference between m^o 
and m n o. We note that the imaginary part of sin 26 is negligibly small, since r^om^o+r^om^o <C 
5m 2 . So, the normalization parameter N of the mixing states is nearly unity. Combining the 
constraints from B(T(3S) — > 7^5) and B(n 10 — > e + e~), we get 

\X d \ = 0.17 ±0.01, m A o~m w o. (33) 

This confirms the results of our straightforward calculation from Eq. (T26]) , but gives a somewhat 
stronger constraint on \X^\. With this constraint, we get 

r(A?) = (9.8 ± 1.1) x 10" 13 MeV, (34) 

which is also in agreement with Eq. (1271) . Furthermore, we get | sin^l 2 = 0.31 ± 0.19. 

It is well known that the decay width of tt° — > 77 agrees perfectly with the SM prediction, 
so it is doubtful that that ir° — > 77 would be compatible with Higgs with a degenerate mass 
m n o. Using the fitted result | sin 6*| 2 = 0.31 ± 0.19 and 

|M(7r /0 -+77)| 2 = i^(|cos^| 2 |^(7r ^77)| 2 + |sin^| 2 |^K^7 7 )| 2 ), (35) 

one can easily observe that 

|^(A°^ 77 )| 2 ~ |^(7r°^77)| 2 (36) 

is needed to give r(7r' — > 77) ~ r(7r° — > 77). However, it would require a too large value of 
\X d \ ~ 10 3 ; therefore, the degenerate case is excluded. 

4 Conclusion 

We have studied the decay 7r° — > e + e~ in the NMSSM and shown that it is sensitive to the 
light CP-odd Higgs boson A\ predicted in the model. The possible discrepancy between the 
KTeV Collaboration measurement P3] and the theoretical prediction of B(tt° — > e + e~) could 
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be resolved in NMSSM by the effects of A\ at the tree level. However, it excludes a large 
fraction of the parameter space of m A o — \Xd\- To further constrain the parameter space, we 
have included bounds from muon g — 2 and the recent searches for A\ from radiative T decays 
performed at CLEO [20] and BaBar [2T]. Combining all these constraints, we have found that 

• B(tt° — >e + e~) and £>(Y — > 'yA®) put strong constraints on the NMSSM parameter Xd and 
m A o. Due to their different dependences on the two parameters, the interesting scenario 
where m A o = 214.3 MeV is excluded, which would invalidate the A\ hypothesis for the 
three HyperCP events |25j. 

• Although these constraints point to a pseudoscalar with m A o ~ m^o and = 0.10±0.08 
(0.17±0.01, 7T° — A\ mixing included) in the NMSSM, such an m A o is excluded by tt° — > 77 
decay. 

In this Letter, we have worked in the limit of Xd 3> X u , i.e., the large tan/3 limit. If 
we relax the limit and take Eqj5] as a general parameterization of the couplings between a 
pseudoscalar and fermions, the u— u— A® coupling should be included. However, its contribution 
is deconstruct ive to the contributions from Xd-, since the 7r° flavor structure is (uu — dd). To 
give a result in agreement with the KTeV Collaboration measurement [16], X u 3> Xd would be 
needed, which would imply possible large effects in ^/(l^) radiative decays. Detailed discussion 
of this issue would be beyond the main scope of our present study. In summary, we could not 
find a region of parameter space of NMSSM with m A o < 7.8GeV in the large tan (3 limit that 
is consistent with the experimental constraints. The HyperCP 214.3MeV resonance and the 
possible 3.3a discrepancy in tt° — > e + e~ decay are still unsolved. Finally, further theoretical 
investigation is also needed to confirm the discrepancy between the KTeV measurements and 
SM predications of ir° — > e + e~ decay. If the discrepancy still persists, it would be an important 
testing ground for NP scenarios with a light pseudoscalar boson. 

The work is supported by the National Science Foundation under contract Nos. 10675039 
and 10735080. 
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